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Multiple Burn Fuel-Optimal Orbit Transfers 
Numerical Trajectory Computation and 
Neighboring Optimal Feedback Guidance 


C.-H. Chuang, Troy D. Goodson, and Laura A. Ledsinger 
School of Aerospace Engineering 
Georgia Institute of Technology, Atlanta, Georgia 30332 

Abstract 

This report describes current work in the numerical computation of multiple bum, 
fuel-optimal orbit transfers and presents an analysis of the second variation for extremal 
multiple bum orbital transfers as well as a discussion of a guidance scheme which may be 
implemented for such transfers. The discussion of numerical computation focuses on the 
use of multivariate interpolation to aid the computation in the numerical optimization. 

The second variation analysis includes the development of the conditions for the 
examination of both fixed and free final time transfers. Evaluations for fixed final time 
are presented for extremal one, two, and three bum solutions of the first variation. The 
free final time problem is considered for an extremal two bum solution, hi addition, 
corresponding changes of the second variation formulation over thrust arcs and coast arcs 
are included. The guidance scheme discussed is an implicit scheme which implements a 
neighboring optimal feedback guidance strategy to calculate both thrust direction and 
thrust on-off times. 

L INTRODUCTION 

The necessary conditions which result from analyzing the first variation of a cost 
functional are widely used. These are commonly referred to as the Euler-Lagrange 
equations. Many problems require additional considerations, for example, the problem 
considered herein, fuel-optimal orbit transfer, requires consideration of Pontryagin's 
Maximum Principle. 
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Many researchers have used the first variation to compute extremal solutions to 
the fuel-optimal orbit transfer problem. Some have used them to apply two-point 
boundary value problem solvers to optimization problems, forming indirect methods. 1>2 > 3 
Others have used a partial set of the conditions to form hybrid indirect/direct methods 
where certain highly sensitive parameters are optimized directly . 4 * 5 However, to the 
knowledge of the authors, few, if any, have made use of the conditions related to the 
second variation of the cost functional. These provide sufficient conditions which, when 
met, declare an extremal solution as a local, weak optimal solution. 

Once the second variation of the cost functional is verified so that it is known 
whether the sufficient conditions are met, the information obtained can then be used to 
implement a guidance scheme. Guidance is defined to be the determination of a way to 
follow an optimal trajectory when presented with obstacles such as environmental 
disturbances or uncertainties in navigation data. Two different types of guidance 
schemes exist: implicit and explicit. Implicit guidance systems are characterized by the 
fact that the vehicle’s motion must be precomputed on the ground and then compared to 
the actual motion. The equations which need to be solved are based upon the difference 
between these measured and precomputed values. The solutions to these equations are 
used in the vehicle’s steering and velocity control. Explicit guidance systems are 
generalized by the fact that the vehicle’s equations of motion are modeled and solved for 
by on-board computers during its motion. The solutions for the equations are solved 
continuously and are used to determine the difference between the vehicle’s current 
motion and its destination. Commands are then generated to alleviate the anticipated 
error. 

Existing guidance schemes have been presented in various papers. An iterative 
guidance scheme which is implemented using a linear tangent steering law is presented 
by Smith . 6 This guidance scheme has been used for the Saturn V and is in currently used 
by the Space Shuttle, the Atlas-Centaur, and the Titan-Centaur. In a paper by Lu 7 * a 
general nonlinear guidance law is developed using two different strategies. One strategy 
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uses optimal control theory to generate a new optimal trajectory onboard from the start, 
while the other uses flight-path-restoring-guidance to bring the trajectory back to the 
nominal. A guidance scheme that is developed using inverse methods for unthrusted, lift- 
modulated vehicles along an optimal space curve is presented by Hough . 8 Linearized 
guidance laws applicable to many different types of space missions are presented by 
Tempelman . 9 These guidance laws are based on fixed and free final time arrivals. 
Naidu 10 presents a guidance scheme applicable to aeroassisted orbital transfers. This 
scheme is developed by implementing neighboring optimal guidance and linear quadratic 
regulator theory. Some interesting techniques for making the neighboring optimal 
guidance converge about the nominal path are introduced in a paper by Powers . 11 

The guidance scheme proposed in this report is an implicit one which implements 
neighboring optimal feedback guidance. An implicit guidance system was chosen due to 
the fact that that type of guidance system handles disturbances well . 10 The neighboring 
optimal feedback guidance was chosen because it inherently uses the nominal solutions. 
Also, it has the advantage of being a feedback system, as a opposed to open-loop 
guidance. 

In this scheme, the initial orbit exit point is assumed to be perturbed from the 
nominal point but the boundary condition, specifying the final orbit, is assumed 
unchanged. The goal is to use the controller to bring the trajectory back to the nominal 
path at some point by using minimal fuel. 

n. First Variation Considerations 

Within this report results are restricted to the planar case, no plane changes are 
considered at this stage of development. The solutions examined in this report satisfy the 
conditions related to the first variation. In the next section, the conditions sufficient for 
declaring a minimizing solution will be checked for these transfers. Some of these 
transfers are multi-bum transfers and in order to simplify initial analysis new nominal 
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solutions have been constructed from these which are single-bum transfers, i.e. the thrust 
is kept on for all time between the initial orbit exit point and the final orbit entry point. 

Only that part of the original trajectory which is contained in the last bum is taken 
to constitute the new extremal solution. This new extremal solution has a fixed initial 
point and fixed transfer time but the final point is only constrained in that it must lie on 
the final orbit. 


n.l. Conditions from the First Variation 

The first order conditions for this problem have been stated many times 12 and will 
be given here only briefly. The optimization problem consists of a cost functional (Eq. 
[1]), state dynamics (Eqs. [2-4]), fixed initial point conditions (Eq. [5]), and boundary 
conditions on the terminal point (Eq. [6]); each of these is expressed below. 
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where r-[x y] T is the radius vector, v=[u v]T is its time derivative, ej is the thrust 
direction, a unit vector, T is thrust magnitude (limited between zero and some maximum 
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value Tjrua), /t is the gravitational constant, g D is the gravitational acceleration at sea- 
level, lap is the specific impulse of the motor. The quantity g 0 Isp is often referred to as 
the exit velocity of the motor. The Euler Lagrange conditions are then 
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where X, y ] T and A. v =CXu X V ] T . The natural boundary conditions are 
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The conditions resulting from applying Pontryagin’s Minimum Principle are 
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where 


H s <0, T = T, 
H s >0, T = 0 


( 16 ) 


H s = 


« So 1 *) 


(17) 


Note that when the derivatives of Hs are zero, singular arc solutions may exist This has 
been checked numerically. 12 

Finally, the free final time problem will also be considered here. For these 
extremal solutions the final, or transfer, time is selected such that the transversality 
condition is satisfied, Le. the Hamiltonian is zero at t=tf. 

tf(,,) = X/r + X. T vL., ( =0 (18 


In a previous paper 12 this problem was given as a maximization problem. To conform to 
the convention used for the second variation 13 , it is now stated as a minimization 
problem. If an extremal solution to the maximization problem is given as state time 
history x(t), Lagrange-multiplier time history X(t), and Lagrange multipliers v, associated 
with boundary conditions, then the extremal solution of the above minimization problem 
isx(t), (-l)*A.(tj,and(-l)*v. 

Additionally, it makes more sense in the guidance problem to consider the control 
as the angle 6, rather than individual components of a unit vector. This simplifies 
analysis because the control is now a scalar. Equation [7] must now be restated as 

tan(d) = ( 19 ) 
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TL2. Extremal Solutions 

All quantities associated with the solutions presented here have been 
nondimensionalized so that fx = 1 and will be presented here in that form. The relations 
used to nondimensionalize are given below. 
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where r* andn* are indicated in the tables for each case of the extremal solutions. 

Each of the transfers given below have both the initial orbit exit point and final 
orbit entry points free. However, for the guidance problem it makes more physical sense 
to consider the initial orbit exit point as fixed and equal to the optimal choice, for it 
cannot be updated once the transfer has begun. 

The last bum of any multi-bum transfer below may also be taken as a complete 
transfer unto itself. The initial point can be chosen as the one at the very first instant (or 
shortly thereafter) of thrusting for the last bum. The final orbit exit point must remain 
unchanged. Obviously, for these choices the natural boundary conditions for final orbit 
entry point are still satisfied. This new transfer can then be considered as a new extremal 
solution, though to an orbit transfer problem with a different initial orbit. 

Figure 1 shows a one bum ascent extremal solution. This trajectory is a transfer 
leaving an orbit with a semimajor axis <2=1.069, eccentricity e=0.02633, and argument of 
perigee cd= -50°. The transfer ends at a nearby orbit with a=1.038 and e=0. Other 
pertinent transfer data are given in Table I. This transfer was produced by shortening the 
time of a two-bum transfer until the coast arc between them vanished. This transfer is 
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therefore both a minimum mass and minimum time extremal solution because mass and 
time have an affine relationship in the one-bum case. 

Figure 2 displays a two-bum transfer, in fact a descending transfer, from an orbit 
with elements a=3.847, e =0.02378 to a final orbit with elements a=1.5, e=0.3333. The 
apses of the terminal orbits are aligned and lie on the 'X' axis of the figure. The initial 
mass is 1.608 and the final mass is 1.1547. Other pertinent transfer data are given in 
Table EL This transfer has the transversality condition converged, therefore it is a 
candidate fuel-optimal free final time solution. By the same right, it can also be 
considered a candidate optimal solution for the fixed final time problem. 



- 0.2 0 0.2 0.4 0.6 0.8 1 1.2 

X 

Figure 1: One-Burn Extremal with Fixed Final Time. 
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Table II Parameters of Transfer in Figure 2. 


A three-bum transfer is shown in Figure 3. Since this transfer is between orbits of 
increasing semimajor axis, it will be referred to as an ascending transfer. The initial orbit 
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has elements a= 2.239, e=0.1 160, and gj=- 85.94°. The elements of the final orbit are «=7, 
e=0.7332, and w=l 14.6°. Other pertinent transfer data are given in Table HI 
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Figure 3: Three-Burn Extremal Orbit Transfer with Fixed Final Time . 


BUB 

0.3898 

ai= 


El 

-85.94° 

Of= 


E 

114.6° 

E2H 

0.01386 

vm 

gheb 

F9 

85.00 

ef= 

MJLJiM 

E 

0.6056 

BS 

6378 km 

ea 

EHH 

| 


Table III. Parameters of the transfer shown in Figure 3. 
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EL Checking the Second Variation 

Extensive derivation of the conditions for the second variation of the cost 
functional has already been detailed in Ref [13]. Equations given in this section are the 
results of applying this work to our problem. 

Considering the second variation of the augmented cost functional, /, a new 
optimal control problem can be stated. In this new problem, the state is Sx , the control 
5u, and the Lagrange-multipliers are SX and dv.. Thus the new cost function is 


S 2 7 = 
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Sx = fJSx + fJSu 
Sxta^Sxo 


(26) 
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where x=[r T v 1 m] T and u=0. 

In general, neighboring optimal feedback guidance allows the designer to consider 
changes in final boundary conditions. We consider no such changes, assuming that the 
destination orbit was accurately planned well in advance. Formulation will be made 
below for both the fixed and free final time cases. 

Evaluating the terms in Eq. [26-28], for orbital transfer, the partial derivatives for 
the dynamics and for the Hamiltonian are: 
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The fixed and free final time problems have the following equations in common.* 


Sx = A(t)Sx-B(t)a 
<5X = -C(t)5x - A T (t)5X. 


(34) 

(35) 


where 


A(t) = f,-f B H;:H lu 
B(t) = f 0 H^f B 
C(t) = H„-H xu H;X 


(36) 

(37) 

(38) 
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For a multiple-bum solution, one finds that Hoe becomes zero during coast arcs. 
This makes it impossible to solve for the change in control, 56. However, since the thrust 
is off during a coast arc, it physically makes no difference what choice is made for the 
control. Therefore, 56 may be chosen as zero and simpler expressions for A(t), B(t), and 
C(t) can be written as 


A(t) = f, 

(39) 

B(t) = 0 

(40) 

C(t) = H„ 

(41) 

Using the sweepback method for nonlinear terminal constraints, as is the case for this 
development, the form for 8X and 5\jr are assumed as 

<5X(t) = P(t)5x(t) + S(t)<fv 

(42) 

5\\f = S T (t)5x(t) + V(t)dv 

(43) 

which allows the solution for dv to be written as 


rfv = V -1 (t 0 )[5v-S r (t 0 )5x(t 0 )] 

(44) 

As mentioned above, 5\|/=0 will be considered here. The boundary condition equations 
are given by: 

P(t f ) = [4»„+(v T V x ) i ] t _ tf 

(48) 
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where in the development for the orbital transfer these are: 
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Following from the assumptions expressed as Eqs. [46-47], the following nonlinear 
equations for P, S, and V must be integrated backwards. The results will be used to 
check the sufficient conditions governing a minimizing solution. 


p=-pa-a t p+pbp-c 


(54) 


s = -(a t -pb)s 


(55) 


v = s t bs 


(56) 


To satisfy the sufficient conditions, Hee > P, S, and V must be such that 


convexity condition: H^t) > 0 for t 0 < t < t f 


(57) 


normality condition: V~ ] (t) exists for t 0 < t < t f 


(58) 


conjugate point condition: P(t) - S(t)V _1 (t)S T (t) finite for t c £ t < t f 


( 59 ) 


The convexity condition is satisfied for any transfer satisfying the choice of 
control specified by the Euler-Lagrange equations. This can be seen by noting that 
Eq.[32] is positive definite, irrespective of the time history for the Lagrange multipliers. 
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m.l. Numerical Results for Fixed Final Time 

The results discussed in this section were obtained for nominal solutions with 
fixed transfer time. The eigenvalues of the V(t) matrix are plotted in Figure 4 for the 3 - 
bum extremal solution. Note that V(t) is not negative definite, one of the eigenvalues is 
zero for all time and the other two eigenvalues are positive. Therefore, the normality 
condition is violated. Furthermore, the conjugate point condition cannot be checked. It 
must be concluded that the 3-bum extremal does not satisfy the sufficient conditions and 
cannot be considered an optimal solution for fixed final time. Figure 5 shows the 
eigenvalues of V(t) over the one-bum extremal constructed from this solution in the 
manner described earlier. The same conclusions must be made for this transfer. 
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Figure 5: Plot of Eigenvalues of V(t) for Last Burn of Three-Burn Extremal. 

The eigenvalues of V (t) for the single-bum transfer are shown in Figure 6. This 
plot is again made for the two-bum extremal in Figure 7 and a one-bum constructed from 
it in Figure 8. These figures all show similar results, namely that V(t) is not negative 
definite, but positive semidefinite. The situation has been repeated, namely that the 
normality condition has been violated and the conjugate point condition cannot be 
checked. Therefore, none of the extremal solutions with fixed final time given in this 
report satisfy the sufficient conditions for a minimizing solution. 
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time (nondimensional) 


Figure 6: Plot of Eigenvalues for One-Burn Extremal. 



time (nondimensional) 

Figure 7: Plot of Eigenvalues of V(t) for Two-Burn Extremal. 
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Figure 8: Plot of Eigenvalues of V(t) for Last Burn of Two-Burn Extremal. 


m.2. N UMERICAL RESULTS FOR FREE FINAL TIME 

When the final time is unspecified, a new condition, the transversality condition, 
must be satisfied by the nominal solution. This condition is expressed in Eq. [60a] 

o(x, u , v 4_ i( = (^r +L ) t =o 

where <I> = <|>(x,t) + v T y(x,t) 
d<h _ . 

~di~lk + lk X 

This slightly complicates the process of checking the sufficient conditions. The 
sweepback method can be used with some additions. Three differential equations and 
thus three boundary conditions must be added to those for P, S, and V. 


(60a) 

(60b) 

(60c) 


m = -(A T -PB)in, 



(61) 
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n = S T Bm , 


a - m T Bm, 



These additions are used to form P , S , and V matrices as follows: 


P = P 


S = S- 


mm 

a 

mn r 
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V = V- 


nn 
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(62) 

(63) 


(64) 

(65) 

( 66 ) 


The equations for dv and 5 A. change by substituting P, S, and V for P, S, and V 
respectively, giving 


A. = V-'(t„)[5 V -S'(t„)«x(g] 

(67) 

«.(t) = P(c)&(t)+5(t)<fv 

(68) 


Note again, however that S\\j has been assumed zero. Now, the sufficient conditions 
based on the second variation with free final time are: 


convexity condition: (t) > 0 for t 0 < t < t f 


(69) 


V _1 (t) exists for t 0 < t < t. (70a) 

normality condition: 

a -1 (t) exists for t 0 <t <t^ (71b) 


conjugate point condition: P(t) - S(t) V"‘(t)S T (t) finite for t„ < t < t f (72) 
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The eigenvalues of V are plotted in Figure 9. Figures 10-12 plot the elements of 
the conjugate point condition matrix. Figure 13 is a plot of a(t). Figure 9 shows that V 
is positive definite in the required interval. Figure 13 shows that ce(t) is negative definite 
in the required interval. Since the normality condition requires that the inverse of V and 
a(t) exists in the interval, this solution is normal. Figures 10-12 show that the conjugate 
point condition is satisfied. The elements are bounded in the required interval and grow 
asymptotically at the final time. Therefore, this solution satisfies the sufficient conditions 
for minimizing the cost functional with free transfer time. 



Figure 9 Plot of Eigenvalues of V (t) for Two Burn Extremal. 
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Figure 10: Plot of Elements of Conjugate Point Condition Matrix for Two Burn 

Extremal. 



Figure 11: Plot of Elements of Conjugate Point Condition Matrix for Two Burn 

Extremal. 
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Figure 12: Plot of Elements of Conjugate Point Condition Matrix for Two Burn 

Extremal. 


Plot Of cc( t ) 



Figure 13: Plot of a(t) for Two Burn Extremal. 
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Numerical results were also obtained for a one bum case using the same free final 
time solution. Figs. 14 & 15 show that the normality condition is indeed met in that, 
respectively, a(f) exists and V 1 exists 



Figure 14: Plot of a(t) for One Burn Extremal. 
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As seen above in Figs. 10-12, the conjugate point condition is met for the two 
bum case in that the elements of the matrix required for that condition are finite; thus the 
conditions can be met similarly for the one bum case for the same solution. 

IV. Neighboring Optimal feedback Guidance 

Conveniently, construction of a neighboring optimal feedback guidance law uses 
the same information as that required to check the second variation of the cost functional. 
As a result, much of the derivation required of guidance law has been stated already. The 
remaining discussion will describe how to form the feedback control law and adjust the 
characteristics of the bang-bang control in a feedback law. 

The control, 89, for the fixed final time problem can be found using 

5«(t) = JP)Sx + f JSrfv] (7J) 

=-H;;[fI(p-sv-'s r )]«x 


Note that this continuous feedback law has been constructed by estimating dv at each 
instant of time. The feedback law depends on P, S, and V as functions of time. A 
particular advantage of the sweepback method is the solution of P(to), S(t 0 ), and V(to), 
allowing the guidance law to store these values and propagate them forward to the current 
time to calculate the current feedback gain. Propagation of the feedback gain may be by 
integration or more practically by interpolation between stored values. Use of this control 
should keep the trajectory on a neighboring optimal solution and deliver the spacecraft to 
the required orbit in the specified transfer time. 

If the transfer time is not fixed, and was chosen optimally for the nominal 

trajectory. Then the formulation for free final time as stated earlier can be used to obtain 
the feedback guidance law 
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time to calculate the current feedback gain. Propagation of the feedback gain may be by 
integration or more practically by interpolation between stored values. Use of this control 
should keep the trajectory on a neighboring optimal solution and deliver the spacecraft to 
the required orbit in the specified transfer time. 

If the transfer time is not fixed, and was chosen optimally for the nominal 
trajectory. Then the formulation for free final time as stated earlier can be used to obtain 
the feedback guidance law 
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<50(0 = -H-[(flP)5x + fJSJv] 
= -H^[fI(-SV- , S T )j5x 


( 74 ) 


and the change in the final time, dt f , is: 


dt, 



(75) 


Evaluating dtf determines when the thrust will be turned off to complete the transfer. 

The block diagram for the feedback controller needed for neighboring optimal 
feedback guidance is shown in Figure 16. 
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Figure 16: Diagram of Neighboring Optimal Feedback Controller 
Implementation. 13 


In Figure 1 6 A ] (t) is the feedback gain for the 80 equation. 

To determine when the new switching times should be, a variation of the 
switching function must be taken 


dH T = H Tx d\ + H TX dk + H Tg d0 = 0 


(76) 
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Therefore, the equation to find the change in the switching time is 



_ X T S\ - x T 5A 
-2x t A 

In order to implement changes in the switching times it will be necessary to predict future 
errors in the state. The state transition matrix should be sufficient in this matter. Such 
predictions will provide the foresight to make switching times earlier or later when 
necessary. 

V. S imulation Results For The Free Final Time Solution U sing The One Burn 
Case 

The controller was implemented by simulating the one bum case for the free final 
time solution. The simulation corresponds to a forward integration of the states, costates, 
and the assumed variables, P, S, V, m, n, and a from the initial time to the final time. 

A comparison of the nominal and actual trajectories is shown in Fig. 17. (The 
actual trajectory being that generated from the simulation results.) Fig. 18 shows a plot 
of the actual and approximated errors in the trajectories when each state is perturbed from 
a value of 10"4 (actual refers to the difference between the nominal and actual trajectories 
and approximated refers to the integrated error). It is seen that the actual and 
approximated errors are concurrent with one another; however, they do not approach zero 
which implies that stability in this case is not guaranteed. Figure 19 shows the 
approximated error during the backward integration in which the error at the final time as 
set to a small number. This plot seems to show a stable response. In order to examine 
response on a more general basis, the 2-norm of the system transition matrix was plotted 
in Fig. 20. Obviously, if the 2-norm went to zero, response to initial conditions would be 
stable. By this plot, it would seem that in general the response will not be decreasing. 
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Figure 17 



: Plot of the Nominal and Actual Trajectories. 
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Figure 18: Plot of the Actual and Approximated Errors. 
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Figure 19: Plot of the Approximated Errors for Backward Integration. 



Figure 20: Plot of the Norm of the Transition Matrix. 


VL Patched Transfer Method Progress 

Recent work in this research project has been directed toward developing a 
numerical computation scheme that performs well for a wide range of acceleration levels 
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and target parameters. Current effort is in assembling a method referred to here as the 
Patched Method. 

The Patched Method is to consist of two phases: a transfer orbit optimization 
phase and a orbit transfer solver phase. The transfer orbit optimization phase is not so 
much concerned with what values of the Lagrange multipliers are required to take the 
craft from orbit A to orbit B as it is with how much time or fuel is required. The orbit 
transfer solver phase, however, is more concerned with obtaining an accurate 
representation of the transfer and is equally concerned with the values of the Lagrange 
multipliers as it is with the transfer time. Finally, it also seems reasonable to desire a 
method that will search for the optimal solution satisfying the target parameters but will 
also, if that fails, be able to return a sub-optimal solution satisfying the target parameters. 
In other words, it would be better to calculate a sub-optimal solution than obtain no 
solution at all. 

Obviously, the key algorithm is one that can quickly determine the minimum 
transfer time (and fuel requirement) between two given orbits in a single bum. One 
approach that may give satisfactory results is an application of multivariate interpolation. 
Interpolation requires calculation and storage of data ahead of time. Therefore, the first 
question is what needs to be stored? To completely specify a problem the following 
twelve (12) values are required: semimajor axis a(to) and a(t f ), eccentricity e(t Q ) and e(t f ), 
true anomaly v<t 0 ) and v(t f ), argument of perigee cu(t 0 ) and G)(tf), mass m(to), thrust T, 
specific impulse l$p, and transfer time, tf. To specify the problem’s solution storage of 
the Lagrange multipliers A^), Xy(t G ), A u (to), A^U^to) is required. 

The first of the nondimensionalization equations, Eqn. [23], says that a(to) can be 
set to unity (a(t Q )=l) for any orbit transfer problem. A simple choice of coordinate axis, 
aligning it with perigee, will set the initial argument of perigee to be set to zero degrees, 
which also works for any orbit transfer problem. Neither eccentricity nor true anomaly 
have such favorable scaling qualities. Additionally, now that the initial values have been 
scaled, the final values cam. , t. The influence of specific impulse can be removed by 
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assuming a constant mass. The mass may be conrected at the end of the bum, so that 

calculation of the following bum is more accurate. Assuming constant mass also 
removes the need to store /U,. 

The influence of the thrust level may be removed by a somewhat restrictive 
assumption, (fir/p) 2 « 1, where 8r is distance between the actual position of the craft and 
a point on a reference orbit with current radius p. This assumption is consistent both with 
low thrust levels, which stay close to a reference orbit for several revolutions, and 
medium thrust levels, which may only stay close to a reference orbit for a few 

revolutions. The advantage is that the assumption linearizes the dynamics and allows the 
solution to be written as 
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Where 0(t,t 0 ) is the state transition matrix for the homogeneous solution. If the initial 
conditions are set to zero, then the solution is linear with respect to the thrust. Now, one 
solution can easily be scaled for any thrust level; however, the resulting solution must 
satisfy the assumption. The orbital elements can then be determined using a Jacobian 
matrix, as shown in Eq. [79] and easily obtained by taking partial derivatives of equations 
used to convert Cartesian coordinates to orbital elements. 
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The number of parameters required to specify a given transfer are reduced to 
seven (7). <5a(tf), e(t 0 ), <5e(tf), &u(tf), v(to), 5v(tf), and if the transfer time stored with this 
data is the minimum time required for the transfer. To store the solution, it is still 
required to know all the Lagrange multipliers. 
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Before proceeding much further with this discussion a few words should be said 
about multivariate interpolation. It is assumed the reader is familiar with univariate, or 
single variable, interpolation in which there is only one independent variable and one 
scalar function of that variable, though any number of such dependent variables is 
allowed. Bivariate interpolation is then interpolation involving two independent variables 
and at least one function of these variables. Bivariate interpolation, and this applies 
equally well to multivariate interpolation, is most easily implemented when the values for 
independent variables are evenly spaced in a grid 14 . 

However, in the case of orbital transfer it would be quite difficult to obtain data 
with the orbital elements of the target orbit evenly spaced because this would require an 
iterative solver for each data point. If that process were easy, there would be no need for 
this approach. On the other hand, it is relatively easy to obtain a grid with the Lagrange 
multipliers, orbital elements of the initial orbit, and the transfer time evenly spaced. The 
equations of motion can then be integrated and the orbital elements of the destination 
orbit are known. The difficulty associated with the unevenly spaced grid is evident when 
one has values for the elements of the target orbit and wants to obtain the values of the 
Lagrange multipliers and the required transfer time. Currently, both types of grids are 
being considered. 

The spacing of the grid is another issue altogether. The spacing of the grid, or its 
density, will determine the accuracy of the estimated minimum transfer time. Since 
speed of the algorithm and storage space required for the software are always important, 
the grid will need a somewhat wide spacing. There are 7 values to store for each point in 
the grid; whichever of the other 4 variables are evenly spaced, it is easy to determine their 
values though proper indexing of the grid points. If each variable is allowed n different 
values and 8 bytes are used to store each number in the computer, then the grid will 

occupy 56/t 7 bytes of memory. For the grid to need under one megabyte of space, then 

\ 

n = 4 is required. For n= 5, the grid would occupy just over 4 megabytes of space. More 
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than likely, different spacings of each variable would be most efficient but will this will 
not change the fact that the grid cannot be dense. 

This interpolation will most surely produce a quick, though rough estimate of the 
transfer time between two chosen orbits. The next phase of the method is to obtain 
accurate solutions for the transfers between these orbits. Estimates for the Lagrange 
multipliers can be also obtained through the interpolation. These can then be used as an 
initial guess for a numerical solver. And, if that fails, a homotopy algorithm can be 
initiated from a nearby grid data point since that data point is already an accurate 
solution. 

vn. Conclusions 

Concerning the calculation of optimal transfers, the current direction has been 
elaborated upon, which is to test numerical methods within the framework of the Patched 
Method. Some work from previous reports and borrowing techniques from the literature 
will be incorporated along with the discussed methods. Results from this work are 
forthcoming. 

A few conclusions, which lend themselves to study in current research work, can 
be made from the analyses presented here. If there are no algorithm mistakes, then it can 
be concluded that the extremal solutions examined may not be locally optimal solutions 
for fixed transfer time. However, some considerations necessary for an accurate 
examination of the second variation may have been overlooked. Ongoing research work 
is in examining why these conditions are not met. 

It was found that the sufficient conditions were satisfied for the free final time 
problem. Software has been developed in order to simulate neighboring optimal 
feedback guidance. Currently, this software is not producing stable solutions. The issue 
of stability of the response must be investigated further and is an area of current research 
endeavors. 
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